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NOMENCLATURE 

Symbol :      Definition: 

A  Closed   loop   dynamics   matrix   for   the 

linearized  system 

a  Control  surface  coordination  gain 

b(x)  Local  beam  of  the  hull 

B  Vehicle  buoyancy 

B  Control  matrix  in  state  space 

CD  Quadratic  drag  coefficient 

clfc2  Coupled  heave  and  pitch  terms 

dq,dw  Cross  flow  drag  terms 

5b  Bow  plane  deflection 

58,5  Stern  plane  deflection 

Iy  Vehicle  mass  moment  of  inertia 

IH,Ip  Cross  flow  drag  terms 

k1,k2/k3,k4  Controller  gains  in  9,w,q,  and  z 

K0  Feedforward  gain 

m  Vehicle  mass 

M  Pitch  moment 

Ma  Partial  derivative  of  M  w.r.t.  a 

0  Pitch  angle 

q  Pitch  rate 

Tc  Time  constant 

U  Forward  speed 

U0  Nominal  speed 

w  Heave  velocity 

W  Vehicle  weight 

x  State  vector 

xB,zB  Body  fixed  coordinates  of  vehicle 

center  of  buoyancy 

xG,zG  Body   fixed  coordinates   of  vehicle 

center  of  gravity 

z  Deviation  off  the  nominal  depth 

zGB  Vehicle  metacentric  height 

Z  Heave  force 

Z=  Partial  derivative  of  Z  w.r.t.  a 
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I .    INTRODUCTION 

The  fundamental  goal  of  submarine  control  is  to  reach  and 
maintain  ordered  depth.  Experimental  designs  involve  expensive 
model  testing  such  as  Darpa  Suboff  Model  (DTRC  Model  5470) 
[Ref.  6].  Much  research  has  been  done  in  depth  control  of 
submarines  [Ref.  3,5].  Our  goal  is  to  develop  an  analytic 
method  to  determine  the  stability  properties  of  a  design. 

The  stability  of  a  design  will  have  a  significant  impact 
on  its  responsiveness.  A  vehicle  with  a  large  margin  to 
instability  will  not  be  very  responsive.  The  problem  becomes 
one  of  determining  how  close  to  the  margins  we  can  get  without 
a  total  loss  of  stability.  By  expanding  the  scope  of  our 
research  to  include  nonlinear  terms  we  are  able  to  define  the 
limits  of  stability  and  therefore  margins. 

Previous  studies  analyzed  stability  properties  of  the 
system,  specially  static  bifurcations  [Ref.  2]  and 
bifurcations  to  periodic  solutions  [Ref.  1].  The  latter  study 
which  is  used  as  a  basis  for  this  work,  was  restricted  to 
level,  zero  trim  flight  paths. 

The  purpose  of  this  thesis  is  to  develop  a  program  for 
finding  the  limits  of  stability  for  an  out  of  trim  submarine 
at  moderate  and  high  speeds.  These  limits  are  mapped  using  a 
Hopf  bifurcation  analysis  program  included  in  the  Appendix. 


II.    EQUATIONS  OF  MOTION 

The  motion  of  the  submersible  in  the  vertical  plane  can  be 
modeled  by  four  coupled  nonlinear  differential  equations  for 
pitch  rate  (q)  ,  heave  velocity  (w)  ,  pitch  anqle  (0)  and  heave 
(z).  With  a  body  fixed  coordinate  frame  at  the  vehicle's 
geometric  center  ,  we  can  express  Newton's  equations  of  motion 
as 


m  (w  -  Uq  -  ztf  2  -  Xaq)  =  Z fl  +  Z^w  +  Z^66  +  Z,6,  +  Z„w  +  Z qq 


-p[CDb(x)± ^-dx 

2    J  I  w  -  xa  I 


e  -q 


z  =  -  Usin  0  +  w  cos  0 


(2.1) 


I   -mxQ(\v  -Uq)  -  zQwqm   =  -xGBBcos  0  -  zQBBsin  0  +Mft  6f 
Mhbb  +  M.q  +  Mjv  +  Ma  +  M  w  +  -p  f  CD&fx,>  j ^-  xdx 


(2.3) 


(2.4) 


Equations  2 . 1  through  2 . 2  can  be  written  in  a  more  compact 
form  as, 


2\, 

s 


w  =  djjUw  +  a]2Uq  +  a]3zGBsin  6  +  a]3xGBcos  6  +  bjU    6 
b2U2bb+dJw,q)  +  Cl(w,q) 


q  =  a2J  Uw  +  a22U q  +  a23zGBsin  9  +  a23xGgcos  0  +  b}  U    6f 
b2u2(>b+  dq(w>i) +  c2(w>4) 


where , 


Dv=  (m-ZJ(Iy-M4)- (mxG+Z4)(mxG  +  MJ 

anDv  =  (Zw-2CDE0Utan  dQ) (Iy-  M f  *  (mx Q+ZJ  (Mw  +  2C DE ^ tan  6^ 

a12Dv=  (m+Zq+2CDE}Utan  dQ)(Iy-M4)  + 
(mxG+  Z.)(M-  mx G-mz GUtan  QQ-  2CDE2Utan  QQ) 

a21Dv  =  (M^  2CDEjUtcm  QQ)(m  -  ZJ  *  (mx  G  +  M»(Z„-  2CDEQUtan  6^ 

a22Dv  =  (M   -  mxG-  mzGUtmdQ-  2CDE2UtmQQ) 
(mxG  +  M^)(m  +  Z^+  2CDExU\m  60) 

a23Dv=  -(m  +  ZJB 

al3Dv  =  -(mxG+Z4)B 

*>Pv=(Iy-Mt)Zh-(G+ZJMh 


(2.5) 


(2.6 


(2.7) 


b2Dv  =   (m~  ZJM6  +  (mXG+M*>  Z6 

dq(w.q)Dv=  (m-ZJIq+(mxG  +  MJI„ 
dJn.q)Dv  -  (Iy-Mq)Iw-  (mxG+  Zq)Iq 


Cj(w,q)Dv  =  (Iy-Mq)mzGq2-  (mx Q+  Z q)mzQwq 


c2(w.q)Dv  =  -  (m  -  Z^)mzGwq-  (mx  G  +  M^)mzGq  2 


In  these  equations  the  submersible  is  assumed  to  be 
neutrally  buoyant  (W=B) ,  and  statically  stable  ( zG  >  zB)  .  Here 
we  can  assume  zB  to  be  zero,  hence   zGB  =  zG 

At  steady  state  the  cross  flow  drag  integral  terms  IH  and 
Ip  have  the  form, 

IH=  -CDw\w\  Cb(x)dx  IP  =  CDw\w\jb(x)xdx  (2.8) 

From  equation  (2.3)  it  is  seen  that  w  is  equal  to  tan0o  at 
steady  state.  The  Jb(x)dx  term  is  computed  numerically  for  the 
SUBOFF  model  as  E0/  and  fb(x)xdx  term  as  E1#  Therefore,  the 
cross  flow  drag  terms  become, 

IH  =   -CD  w   jwl  E0  Ip  =   CDw   jwl  E1  (2.9) 
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Because  we  have  two  rudders  at  the  bow  and  the  stern,  our 
system  of  equations  is  multi-input.  To  reduce  this  system  into 
a  single  input  system  the  linear  combination  of  the  control 
inputs  will  be  modified  into  the  following  form, 

5   =  ds        ,        5b=  a53  (2.10) 

where  a  is  the  control  surface  coordination  gain.  The  value 
of  a  ranges  from  -1  to  1.  The  selection  of  the  value  of  a  will 
allow  the  planes  to  operate  as  desired  for  the  particular 
maneuvering  condition,  i.e.,  a  =  0  for  no  bow  plane  control, 
a  =  -1  for  bow  plane  and  stern  plane  opposed  to  each  other, 
yielding  the  maximum  pitch  moment,  and  a  =  1  for  bow  and  stern 
plane  control  in  the  same  direction,  yielding  the  maximum 
heave  force. 


.Vcal.zal 
Deoch 


He: izcncdi 


S3 


'J 


FIGURE   2.1   Geometric   representations   of   the   basic 
definitions . 


III.    CONTROL  LAW 
A.   INTRODUCTION 

The  control  design  problem  can  be  expressed  in  state  space 

as  follows, 

(3.1) 
x  =  Ax  +  J36  v    ' 

where  the  state  vector  is 


x  = 


(3.2) 


Equation  3.1  in  our  case  is, 


■    - 

e 

\v 

= 

4 

z 

7  7  7  J 

al3zGB~blu   kl    auu~bjU   k2    cijji-bjU   k3    -bjU   k4 
a23ZGB'b^2kl     a21U-b^2k2    a22U~b2u2k3     ~b/*2k4 


0 


0 


(3.3) 


Our  aim  is  to  find  a  controller  which  will  assure  us  a 
stable  closed  loop  system.  The  only  control  input  is  the  dive 
plane  angle,  5. 


B.   FEEDBACK  CONTROL 
1 .   Pole  Placement 

The  full  state  feedback  controller  is  a  linear  function  of 
the  states  and  has  the  form, 

5  =   -   K  x  (3.4) 

where  K  is  the  vector  of  feedback  gains  which  are  to  be 
determined  in  order  to  give  the  desired  closed  loop  system 
dynamics.  Substituting  equation  2.12  into  2.10  yields, 

x  =  (A  -  BK  )  x  (  3  • 5  ) 

The  feedback  gains  K  must  be  chosen  such  that  A  -  BK   has  the 

desired  eigenvalues.  The  actual  characteristic  equation  of  the 

closed  loop  system  is  given  by, 

det(A   -   BK   -    si)    =   0  (3.6) 

The  required  values  of  K  are  obtained  by  matching  coefficients 

in   the   two  polynomials   of   the  actual   and  the  desired 

characteristic  equations.  Equation  3.5   becomes, 


w 


0           0  10 

a13ZGB  allu  al?      ° 

a2SZGB  a21U  a22U     ° 

-u          1  0       0 


0 

0 

w 

+ 

bju2 

<l 

b2u2 

z 

0 

(3.7) 


The  characteristic  equation  of  the  closed  loop  system  is, 


det 


-  s 


0  JO 

")                                          0  7                                7 

OjfGB~ b jU  kj    djjti-bjU  k2~s  a12u-b}u   k3       -bjU  k4 

a23zGB-b2'i2kl        a2l"-b^2k2  a22U  ~  h  ^  3~  S     ~b2,i  % 

1  0  -s 


-u 


(3.8) 


which  reduces  to, 


s4  +  (A2k2  +  A3k3-E])s3+(-B}k1-  B2k2 -  B3k3 -  B4k4 -  E2)s2  + 
(-C1k]-C2k2-C3k3-E3)s+  (-Dlk4-D2k4)  =  0 


(3.9) 


where, 

A2  =   -B4   =  bx  u2 

A3  =    -B1    =  b2  u2 

B2  =    (a22   b,    -   a12   b2)    u3 

B3  =   Ci    =    (a n   b2   -   a21  bj    u3 

C2  =  Di    =    (a23  bl  -a13  b2  )    zGB   u2  (3.10) 

C4  =    (a22  bx    +  b2   -   a12  b2)    u3 

D2  =    (an   b2   -   a21  b:)    u4 

Ei  =    (au    +   a22)    u 

^2      ~      a23ZGB      +       (ai2a21      ~      311      322)       U 
E3     =      (313     321      "     ail     a23)      ZGB     U 

Now,  let's  assume  that  we  want  to  place  the  closed  loop  poles 
at  -plf  -p2,  -p3,  -p4  to  have  the  desired  system  response.  Then 
the  desired  characteristic  equation  is, 

4  3  2  n  (3.11) 

S      +  CtjS      +  «25   +  a3S   +  a4  =   " 


where, 

ai    =  Pi    +  P2   +  P3   +  P4 

&2  =  P1P2   +  PiPj   +Pi  P4    +  P2Pj   +  P2P4    +  P3P4  (3.12) 

<*>  =  Pi  P2  P3    +    P:  p2  p4    +    Pi  p3  p4    +    p2  p3  p4 

°<4    =  Pi  P2    P3    P4 

The  feedback  gains  can  now  be  computed  by  equating  the 

coefficients  of  equation  3.9  and  3.11, 

A2   k2   +  A3   k3   =   -a2    -   E, 

B2   k2   +  B2   k2   +  B3   k3   +  B4   k4   =  a2   +  E2  (3.13) 

Cj  k2    +  C2   k2   +   C4   k4  =   a3   +  E3 

(D,    +  D2)    k4   =  a4 

We  established  the  method  for  placing  the  poles  of  the 
system,  but  we  also  need  to  know  the  desired  locations  of  the 
poles. 

2.   Pole  Location  Selection 

In  a  typical  second  order  system  control  law  design  , 
transient  response  specifications  are  given.  This  results  in 
an  allowable  region  in  the  s-plane  where  the  desired  location 
of  the  poles  can  be  obtained.  For  higher  order  systems  the 
concept  of  dominant  roots  can  be  employed.  In  selecting  poles 
the  physics  of  the  system  must  be  considered.  If  the  poles  are 
too  negative,  a  small  time  constant  will  result,  and  the 
system  may  not  be  able  to  react  that  fast.  If  we  use  big  gains 
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K,  this  also  means  that  the  control  effort  will  be  large.  In 
practice  there  are  limits  based  on  actuator  size  and 
saturation. 

Considering  the  control  law  design  to  stabilize  the 
submarine  to  a  level  flight  path  at  0  =  0  it  is  required  that 
the  submarine  return  to  the  level  flight,  after  some  small 
disturbances  in  9  or  z,  within  the  time  it  takes  for  the 
vehicle  to  travel  three  ship  lengths.  Since  our  model  is  14 
feet  long  and  its  velocity  is  5  feet/sec,  the  required 
recovery  time  is  about  10  seconds.  This  means  the  time 
constant  is  3  seconds  and  the  closed  loop  poles  should  be 
placed  at  approximately   -0.3. 

After  placing  the  poles  using  equation  3.13  the 
control  law  is  found  to  be, 
5   =  0.9917    0  +  0.8333  w  +  0.6026  q   -  0.0351    z  (3.14) 

C.   FEEDFORWARD  CONTROL 

The  previous  discussion  on  feedback  controller  assures 
closed  loop  stability,  but  it  acts  as  a  regulator  in  other 
words  takes  all  the  states  to  zero  values.  If  we  have  constant 
disturbances  or  we  want  to  track  some  reference  value  other 
than  zero  we  can  not  do  this  with  state  feedback  alone. 
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p^ 


r 


reference  signal 


■*Ko 


ieafcr-A/ard  gain 


sum 


_^        x1  =  Ax+du 
y  =  Cx~Ou 


Linearized  modei  of  plant 


K    W- 


vector  of  feedback  gains 


-fe»  yout 
outf. 


FIGURE  3 . 1  Feedback  and  feedforward  control  application  to 
our  linearized  model 


In  the  case  of  non-zero  set  points  or  constant 
disturbances  we  again  need  to  have  the  exact  same  full  state 
feedback  to  ensure  closed  loop  stability.  But  we  also  need  to 
introduce  an  additional  term  to  our  controller  in  the  form 


u  =  -  K  x 


(3.15) 
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where  k0  is  the  constant  feedforward  term.  k0  is  given  as 
[Ref.  4], 

k0   =  H^fO)    x0  (3.16) 

where  x0  is  the  reference  values  of  the  states  and  Hc-1(s)   is 
the  closed  loop  transfer  function, 

H~l(s)    =  C    (si   -  A   +  BR)'1  B  (3.17) 

Another  way  of  getting  k0  is  looking  at  the  steady  state 
equations  of  motion.  In  steady  state  all  the  time  derivatives 
in  equations  2.1  through  2.4  go  to  zero  and  we  have, 

(3.18) 
(3.19) 


w  =  tan  00 


Zbb+Z»w-IH=  ° 


-(xG-xB)BcosQ0-zGBBsinQ0  +  Mb6+Mvw  +  Ip  =  0  (3.20) 

If  the  equation  3.19  is  multiplied  with  M6  and  equation  3.20 
with  Z6  and  set  equal  to  each  other  and  plug  in  the  equation 
3.18,  we  have  an  equation  depending  only  on  0. 


(ZyvMb-MJZb)tanB0+xGBZbcos90+zGBZ()smQ0  + 


Where  E0  and  E1  are  the  integral  terms  computed  numerically  for 
the  Suboff  model.  The  steady  state  value  of  60  is 
found  from  equation  3.20  by  using  a  Newton-Raphson  method 
in  Bifurl  program  in  the  Appendix.  Then  we  can  get  5  from 
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equation    3.19, 


6 


IH~ZyfW        -CDEgtanQ0\tane0\-Zww  (3.22) 


*6  Z6 


After  getting  5,  we  easily  find  k0  from  equation  3.15, 

k0  =   6  +  Kx  (3.23) 

A  plot  of  the  steady  state  angle  0O  versus  xG  for  different 
values  of  zG  is  shown  in  Figure  3.2. 
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FIGURE  3.2   Steady  state  pitch  angle  0O  as  metacentric 
height  varies 
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IV.    BIFURCATION  ANALYSIS 

A.   STABILITY 

The  nonlinear  equations  of  motion  in  the  dive  plane  2.1 
through  2.4  can  be  expressed  in  a  compact  form  as  follows, 

x=f(x)  f4-1* 

Where  x  is  the  state  variable  vector  x=[8,w,q,z].  We  know 
that  the  equilibrium  points,  x0  of  the  system  are  defined  by, 

f(x0)    =  0  (4.2) 

This  is  a  nonlinear  system  of  algebraic  equations  and  it 
may  have  multiple  solutions  in  x0,  which  means  that  the 
nonlinear  system  may  have  more  than  one  positions  of  static 
equilibrium.  If  we  pick  one  equilibrium,  x0  we  can  establish 
its  stability  properties  by  linearization.  The  linearized 
system  becomes, 

it-Ax  <4-3) 

where  A  is  the  Jacobian  matrix  of  f  (x)  evaluated  at  x0, 


a-  SL 

dx 


*b  (4.4 
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and  the  state  has  been  defined  to  designate  small  deviations 
from  the  equilibrium  x0, 

x  — ►  x-x0  (4.5) 

In  system  dynamics  as  long  as  all  eigenvalues  of  A  have 
negative  real  parts,  we  know  that  the  linear  system  will  be 
stable.  This  means  that  the  equilibrium  x0  will  be  stable  for 
the  nonlinear  system  as  well.  This  is  in  fact  Lyapunov's 
linearization  technique. 

B.   BIFURCATION 

Values  of  the  nonlinear  system  parameters  at  which  the 
qualitative  nature  of  the  system's  motion  changes  are  known  as 
critical  or  bifurcation  values.  The  phenomena  of  bifurcation, 
i.e.,  quantitative  change  of  parameters  leading  to  qualitative 
change  of  system  properties,  is  the  topic  of  bifurcation 
theory.  Euler  buckling  (Pitchfork  bifurcation) ,  limit  cycles 
(Hopf  bifurcation)  are  common  examples  of  bifurcation. 

Classical  definition  of  stability  states,  that  the  real 
part  of  all  the  eigenvalues  of  the  system  must  be  negative. 
Therefore,  our  initial  investigations  into  the  stability  of 
the  SUBOFF  model  was  to  find  those  eigenvalues  whose  real 
parts  cross  the  imaginary  axis.  We  used  the  bifurcation 
analysis  program,  included  in  the  Appendix,  to  calculate  the 
eigenvalues  of  the  system. 


17 


By  linearizing  the  equations  of  motion,  equations  2.1 
through  2.4  ,  the  state  space  equations  of  the  dynamic  system 
can  be  written  in  the  form, 

x=Ax+Bu  (4,5) 


where, 

u   =  -  K  x  (4.6) 

and  K  is  the  vector  of  controller  gains,  as  calculated  by  pole 
placement  in  equations  3.13.  The  eigenvalues  of  the  system  are 
found  by  solving, 

det|A-BK-sI|  =  0  (4.7) 

In  the  bifurcation  program  a  pseudo-root  locus  method  is 
employed  where  the  time  constant,  Tc,  is  fixed.  The  constant 
Tc  fixes  to  placement  of  the  system  poles  at  a  given  nominal 
forward  speed  U0  and  then  the  model  speed,  U,  is  varied 
incrementally  with  the  system  eigenvalues  calculated  at  each 
speed  increment.  When  the  real  part  of  an  eigenvalue  changes 
sign  between  the  limits  of  a  speed  increment  a  bisection 
method  is  employed  to  find  the  speed  where  the  real  part  of 
the  eigenvalue  equals  zero. 

For  each  point  where  the  real  part  of  an  eigenvalue 
crosses  the  imaginary  axis  the  associated  Tc  and  U  are  plotted 
on  a  bifurcation  map.  This  map  delineates  the  regions  of 
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classical  stability  (all  eigenvalues  on  the  left  hand  plane) 
from  the  regions  of  instability.  A  family  of  bifurcation  maps 
were  generated  by  varying  nominal  speed,  U0,  initial 
stability,  zGB,  and  longitudinal  center  of  gravity/buoyancy 
separation,  xGB  of  the  submersible. 

Figure  (4.1)  shows  a  typical  bifurcation  map  with  its  five 
distinct  regions  [Ref.  1].  Region  I  is  the  area  of  classical 
stability.  In  region  II  there  is  one  real  positive  eigenvalue 
which  is  indicative  of  a  pitchfork  bifurcation.  Pitchfork 
bifurcations  of  this  model  were  previously  examined  by  Reidel 
[Ref.  1].  Regions  III,  IV,  and  V  have  at  least  one  pair  of 
complex  conjugate  eigenvalues  with  a  positive  real  part.  This 
would  indicate  that  there  should  be  an  unstable  oscillatory 
behavior  for  the  model. 

C.   RESULTS  AND  DISCUSSION 

The  classical  stability  region  in  bifurcation  maps  lies 
between  pitchfork  and  Hopf  bifurcation  boundaries.  The  limits 
or  parameters  must  be  defined  for  the  system  designer  prior  to 
starting  the  design.  By  maximizing  the  region  of  stability  we 
can  give  the  designer  the  most  leeway  in  his  work.  There  are 
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FIGURE  4.1  A  typical  bifurcation  map  showing  the  five 
distinct  regions 
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two  parameters  that  we  used  to  change  the  bifurcation  maps, 
the  longitudinal  separation  between  center  of  gravity  and 
center  of  buoyancy,  xGB  and  the  initial  stability,  zGB. 

First  we  look  at  the  change  in  xGB.  In  figures  4.2  through 
4.7  we  plotted  bifurcation  curves  for  different  initial 
stabilities  as  xGB  varies.  We  can  see  that  as  xGB  increases 
the  Hopf  bifurcation  branches  Hi  an  H2  move  towards  higher 
speeds  and  time  constants  and  thus  increasing  the  stability 
area.  The  H3  branch  however  remains  constant. 

The  other  important  point  that  we  observed  is  that  the 
system  becomes  unstable  at  nominal  speed  at  higher  time 
constants.  This  is  unexpected  because  we  are  designing  around 
our  nominal  speed.  A  more  careful  examination  in  the  trimmed 
case  shows  that  the  actual  forward  velocity  becomes,  \/u2+w2. 
Therefore  the  system  may  become  stable  at  a  value  of  u  other 
than  nominal . 

The  next  parameter  we  examined  was  the  initial  stability, 
zGB.  Figures  4.8  through  4.14  show  the  effect  of  varying  zGB 
from  .2  to  .4  ft  for  different  xGB  values.  The  H3  branch 
remains  constant  while  the  upper  speed  H2  branch  moves  down 
effectively  decreasing  the  area  of  stability.  The  low  speed 
curve  Hi  moves  upwards  and  increases  the  low  speed  area  of 
stability. 
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FIGURE   4.2   Bifurcation  map   as   xg   changes    between   xg=0    and 
xg=-.2,    zg=.2. 
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FIGURE  4.3  Bifurcation  map  as  xg  changes  between  xg=0  and 
xg=.2,  zg=.2. 
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FIGURE  4.4  Bifurcation  map  as  xg  changes  between  xg=0  and 
xg=- . 3 ,  zg=.3. 
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FIGURE   4.5   Bifurcation  map   as   xg   changes    between   xg-0    and 
xg=.3,    zg=.3. 
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FIGURE  4.6  Bifurcation  map  as  xg  changes  between  xg-0  and 
xg=- . 3 ,  zg=.4. 
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FIGURE  4.7   Bifurcation  map  as   xg  changes   between  xg=0    and 
xg=.3,    zg=.4. 
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FIGURE    4.8       The    effects    of    changing    zg    on    the    bifurcatic 
maps ,    xg= . 3 . 
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FIGURE    4.9       The    effects    of    changing    zg    on    the    bifurcation 
maps ,    xg= . 2 . 
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FIGURE  4.10   The  effects  of  changing  zg  on  the  bifurcatio 
maps ,  xg= . 1 . 
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FIGURE    4.12      The    effects    of    changing    zg   on   the    bifurcatio 
maps,    xg=-.l. 
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V.   CONCLUSIONS  AND  RECOMMENDATIONS 

A.   CONCLUSIONS 

Hopf  bifurcation  analysis  is  a  very  useful  design  tool  in 
the  design  and  evaluation  phase.  Hopf  bifurcation  analysis  and 
an  identification  program  that  can  evaluate  the  hydrodynamic 
coefficients  for  the  submersible  vehicle  will  be  very  useful 
and  save  money  and  time  by  reducing  the  amount  of  model 
testing.  An  effective  set  of  control  system  parameters  can  be 
generated  in  this  process  that  will  be  optimal  for  the  final 
design  of  the  submersible. 

This  type  of  analysis  can  set  the  limits  of  the  ranges  of 
important  parameters  such  as  metacentric  height  and 
longitudinal  separation  of  buoyancy/gravity  centers.  As  we 
have  seen  changes  in  these  two  parameters  can  have  dramatic 
effects  on  stability.  It  was  found  that  the  moderate  speed 
region  of  stability  increases  with  increasing  metacentric 
height.  The  same  is  not  true,  however,  for  high  speeds.  The 
longitudinal  separation  of  center  of  gravity/buoyancy  can  have 
a  profound  effect  on  stability.  It  was  found  that  the  vehicle 
may  be  unstable  even  at  nominal  speed.  This  was  attributed  to 
the  fact  that  at  high  trim  angles,  the  feedback  gains  which 
are  computed  at  zero  trim,  can  no  longer  guarantee  stability. 
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B .   RECOMMEND AT I ONS 

The  bifurcation  analysis  program  should  be  expanded  to 
evaluate  the  performance  of  the  submarine  including  effects  of 
external  forces  such  as  wave  effects,  currents,  and  free 
surface  effects. 
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APPENDIX    -   BIFURCATION  ANALYSIS    PROGRAM 


C      PROGRAM  BIFUR1  FOR 
C      BIFURCATION  ANALYSIS 
C      PARAMETERS  ARE:  TC  VS.  U 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  Kl,K2,K3,K4,L,MQDOT,MWDOT,MQ,MW,MDS,MDB,MD, 

&  MASS,IY,P1,P2,P3,P4,XGB,ZGB 

DIMENSION  A(4,4),FV1(4),IV1(4),ZZZ(4,4),WR(4),WI(4),XL(25), 

&  BR(25),VEC0(25),VEC1(25),VEC2(25) 

COMMON  P1,P2,P3,P4 
C 

OPEN  (1 1,FILE=,BIF1.RES,,STATUS=,NEW,) 

OPEN(12,FILE=,BIF2.RES,,STATUS=,NEW,) 

OPEN(13,FILE='BIF3.RES,,STATUS=rNEW,) 
C  NUMERIC  INFO  OF  DARPA  SUBOFF  MODEL 

WEIGHTS  556.2363 

BUO    =1556.2363 

L      =13.9792 

IY     =561.32 

G      =32.2 

MASS  =WEIGHT/G 

RHO    =1.94 

XB     =0.0 

ZB     =0.0 

CD     =0.4 
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CD     =0.5*CD*RHO 


C 


WRITE  (*,*)  'ENTER  MIN,  MAX,  AND  INCREMENTS  IN  Tc  (nondim)' 

READ  (*,*)  TCMIN,TCMAX,ITC 

WRITE  (*,*)  'ENTER  MIN,  MAX,  AND  INCREMENTS  IN  U  (nondim)' 

READ  (*,*)  UMIN,UMAX,IU 
C       WRITE  (*,*)  'ENTER  NOMINAL  SPEED' 
C       READ  (*,*)U0 

WRITE  (*,*)  'ENTER  XG  AND  ZG' 

READ  (*,*)XG,ZG 

U0=9 
C 

ZGB=ZG-ZB 

XGB=XG-XB 

TCMIN=TCMIN*L/UO 

TCMAX=TCMAX*L/UO 

UMIN  =UMIN*U0 

UMAX  =UMAX*U0 
C    HYDRODYNAMIC  COEFFICIENTS 

ZQDOT=-6.3300E-04*0.5*RHO*L**4 

ZWDOT=-1.4529E-02*0.5*RHO*L**3 

ZQ    =  7.5450E-03*0.5:|tRHO'|eL**3 

ZW    =-1.3910E-02*0.5*RHO*L**2 

ZDS  =-5.6030E-03*0.5*RHO*L**2 

ZDB  =-5.6030E-03*0.5*RHO*L**2 

MQDOT=-8.8000E-04*0.5*RHO*L**5 

MWDOT=-5.6100E-04*0.5*RHO*L**4 
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MQ  =-3.7020E-03*0.5*RHO*L**4 
MW  =  1.0324E-02*0.5*RHO*L**3 
MDS  =-2.4090E-03*0.5*RHO*L**3 
MDB  =  2.4090E-03*0.5*RHO*L**3 


XL(1)=0.0 
XL(2)=0.1 
XL(3)=0.2 
XL(4)=0.3 
XL(5)=0.4 
XL(6)=0.5 
XL(7)=0.6 
XL(8)=0.7 
XL(9)=1.0 


XL(10 
XL(11 
XL(12 
XL(13 
XL(14 
XL(15 
XL(16 
XL(17 
XL(18 
XL(19 
XL(20 
XL(21 


=2.0 

=3.0 

=4.0 

=7.7143 

=10.0 

=15.1429 

=16.0 

=  17.0 

=18.0 

=  19.0 

=20.0 

=20.1 
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XL(22)=20.2 
XL(23)=20.3 
XL(24)=20.4 
XL(25)=20.4167 
DO  102  N=  1,25 

XL(N)  =  (L/20.)*XL(N)  -  L/2. 
102  CONTINUE 
BR(1)=0.0 
BR(2)=0.485 
BR(3)=0.658 
BR(4)=0.778 
BR(5)=0.871 
BR(6)=0.945 
BR(7)=1.010 
BR(8)=  1.060 
BR(9)=1.18 


BR(10 
BR(11 
BR(12 
BR(13 
BR(14 
BR(15 
BR(16 
BR(17 
BR(18 
BR(19 
BR(20 


=1.41 

=1.57 
=1.66 
=  1.67 
=  1.67 
=1.67 
=  1.63 
=  1.37 
=0.919 
=0.448 
=0.195 
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BR(21)=0.188 

BR(22)=0.168 

BR(23)=0.132 

BR(24)=0.053 

BR(25)=0.0 

DO  104  K=  1,25 
VEC0(K)=BR(K) 
VEC 1  (K)=XL(K)*BR(K) 
VEC2(K)=XL(K)*XL(K)*BR(K) 
104  CONTINUE 

CALL  TRAP(25,VECO,XL,EO) 

CALL  TRAP(25,VEC1,XL,E1) 

CALL  TRAP(25,VEC2,XL,E2) 
C 

ALPHA=0.0 

ZD=ZDS+ALPHA*ZDB 

MD=MDS+ALPHA*MDB 

C      CALCULATING  THE  SS  PITCH  ANGLE  0O 
C      WITH  NEWTON  RAPHSON  METHOD 

Pl=  ZW*MD  -  MW*ZD 

P2=  XG*BUO*ZD 

P3=  ZG*BUO*ZD 

P4=  CD*(MD*E0  -  ZD*E1) 

WRITE(*,*)P1,P2,P3,P4 

EPSI=.  00000001 

THETA0=0 
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IF(XGB.GT.0)P4=-1*P4 

18  DO  19  1=1,2000 
FT=FUNC(THETA0) 
DFT=DFUNC(THETAO) 
DELT=FT/DFT 
THETA0=THETA0-DELT 
IF(ABS(DELT)-EPSI)  20,20,19 

19  CONTINUE 
FT=FUNC(THETA0) 

20  WRJTE(V)  THETAO  ,  FT 

c: 

DV=(MASS-ZWDOT)*(IY-MQDOT)-(MASS*XG+ZQDOT)*(MASS!,:XG+MWDOT) 
A11DV=(IY-MQDOT)*(ZW-2*CD*EO*U*TAN(THETAO)) 
&      +(MASS*XG+ZQDOT)*(MW+2*CD*E1*U*TAN(THETAO)) 

A!2DV=(IY-MQDOT)*(MASS+ZQ+2*CD*E1*U*TAN(THETAO))+ 

&  (MASS*XG+ZQDOT) 

&      *(MQ-MASS*XG-MASS*ZG*U*TAN(THETA0)-2*CD*E2*U*TAN(THETA0)) 
Al  3DV=-(MASS*XG+ZQDOT)*WEIGHT 
B1DV=(IY-MQD0T)*ZD+(MASS*XG+ZQD0T)*MD 
A21DV-(MASS-ZWDOT)*(MW+2*CD*El*U*TAN(THETA0)) 

&      +(MASS*XG+MWDOT)*(ZW-2*CD*EO*U*TAN(THETA0)) 

A22DVKMASS-ZWDOT)*(MQ-MASS*XG-MASS*ZG*U*TAN(THETA0)-2*CD*E2*U 
&    *TAN(THETA0))+(MASS*XG+MWDOT)* 
&  (MASS+ZQ+2*CD*E1*U*TAN(THETA0)) 

A23DV=-(MASS-ZWDOT)*WEIGHT 
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c 


B2DV=(MASS-ZWD0T)*MD+(MASS*XG+MWD0T)*ZD 

A11=A11DV/DV 

A12=A12DV/DV 

A13=A13DV/DV 

A21=A21DV/DV 

A22=A22DV/DV 

A23=A23DV/DV 

B1=B1DV/DV 

B2  =B2DV  /DV 

EPS  =l.D-5 
ILMAX=1500 

DO  1  1=1,  ITC 
WRITE  (*,2001)I,ITC 

TC=TCMIN+(I- 1  )*(TCMAX-TCMIN)/(ITC- 1 ) 
POLE=1.0/TC 
ALPHA3=4.0*POLE 
ALPHA2=6.0*POLE**2 
ALPHA1=4.0*POLE**3 
ALPHA0=    POLE**4 

K4=ALPHA0/((B1*A21-B2*A11)*U0**4+(B1*A23-B2*A13)*ZGB*U0**2) 
A2M=B1*U0**2 
A3M=B2*U0**2 
A0M=-(A1 1+A22)*U0-ALPHA3 
B1M=B2*U0**2 
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B2M=(B2*A12-B1*A22)*U0**3 

B3M=(B1*A21-B2*A1 1)*U0**3 

B0M=(A1 1  *  A22-A2 1  *  Al  2)*U0**2-A23  *ZGB-ALPHA2-B  1  *U0*U0*K4 

C1M=(B2*A1 1-B1*A21)*U0**3 

C2M=(B1*A23-B2*A13)*ZGB*U0**2 

C0M=(A13*A21-A23*A11)*ZGB*U0+ALPHA1-(B2+B1*A22-B2*A12)*K4*U0**3 

K2-C1M*B0M*A3M-B1M*C0M*A3M-C1M*B3M*A0M 

K2=K2/(C1M*B2M*A3M-B1M*C2M*A3M-C1M*B3M*A2M) 

K 1  =(C0M-C2M*K2)/C  1 M 

K3=(A0M-A2M*K2)/A3M 


D0  2J=1,IU 
U=UMIN+(  J- 1  )*(UM  AX-UMIN)/(IU- 1 ) 
A(1,1)=0.0D0 


A(l,2 
A(l,3 
A(l,4 
A(2,l 

A(2,2 
A(2,3 
A(2,4 
A(3,l 
A(3,2 
A(3,3 
A(3,4 
A(4,l 
A(4,2 


O.ODO 
=1.0DO 
=O.ODO 

=-Al  3  *(XGB*SIN(THETAO)-ZGB*COS(THETAO))+B  1  *U*U*K1 
=A11*U  +B1*U*U*K2 
=A12*U  +B1*U*U*K3 

B1*U*U*K4 
=-A23*(XGB*SIN(THETA0)-ZGB*COS(THETA0))+B2*U*U*Kl 
=  A21*U  +B2*U*U*K2 
=  A22*U  +B2*U*U*K3 
=  B2*U*U*K4 
=-U 
=  l.ODO 
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c 


A(4,3)=  O.ODO 
A(4,4)=  O.ODO 

CALLRG(4,4,A,WR,WI,0,ZZZ,IV1,FV1,IERR) 
CALL  DSTABL(DEOS,WR,WI,FREQ) 

IF(J.GT.l)GOTO10 

DEOSOO=  DEOS 

UOO  =  U 

LL=  0 

GO  TO  2 
10      DEOSNN  =  DEOS 

UNN  =  U 

PR=  DEOSNN*DEOSOO 

IF  (PR.GT.O.DO)  GO  TO  3 

LL  =  LL+1 

IF  (LL.GT.3)  STOP  1000 

IL  =  0 

UO  =  UOO 

UN=  UNN 

DEOSO  =  DEOSOO 

DEOSN  =  DEOSNN 
6      UL  =  UO 

UR  =  UN 

DEOSL  =  DEOSO 

DEOSR  =  DEOSN 
U  =  (UL+UR)/2.D0 
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ALPHA  =  (DEOSL-DEOSR)/(UL-UR) 

U  =-  (DEOSL-ALPHA*UL)/ALPHA 

A(l,l)  =  O.ODO 

A(l,2)  =  O.ODO 

A(l,3)  =  1.0D0 

A(l,4)  =  O.ODO 

A(2,l)  =  A13*(XGB*SIN(THETAO)-ZGB*COS(THETAO))+Bl*U*U*Kl 

A(2,2)  =  A11*U  +B1*U*U*K2 

A(2,3  )=  A12*U  +B1*U*U*K3 

A(2,4  )=     B1*U*U*K4 

A(3,l)  =  A23*(XGB*SIN(THETAO)-ZGB*COS(THETAO))+B2*U*U*Kl 

A(3,2)=  A21*U  +B2*U*U*K2 

A(3,3)=  A22*U  +B2*U*U*K3 

A(3,4  )  =     B2*U*U*K4 

A(4,l)=-U 

A(4,2)=  l.ODO 

A(4,3)  =  O.ODO 

A(4,4  )=  O.ODO 

CALLRG(4,4,A,WR,WI,0,ZZZ,IV1,FV1,IERR) 
CALL  DSTABL(DEOS,WR,WI,FREQ) 

DEOSM  =  DEOS 

UM  =  U 

PRL  =  DEOSL*DEOSM 

PRR=  DEOSR*DEOSM 

IF  (PRL.GT.O.DO)  GO  TO  5 
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UO  =  UL 

UN  =  UM 

DEOSO  =  DEOSL 

DEOSN  =  DEOSM 

IL  =  IL+1 

IF  (IL.GT.ILMAX)  STOP  3 100 

DIF  =    DABS(UL-UM) 

IF  (DIF.GT.EPS)  GO  TO  6 

U  =  UM 

GO  TO  4 
5      IF  (PRR.GT.0.D0)  STOP  3200 

UO  =  UM 

UN  =  UR 

DEOSO  =  DEOSM 

DEOSN  =  DEOSR 

IL  =  IL+1 

IF  (IL.GT.ILMAX)  STOP  3100 

DIF  =  DABS(UM-UR) 

IF  (DIF.GT.EPS)  GO  TO  6 

U  =  UM 
4      LLL  =  10+LL 

WRITE  (LLL,*)  U/U0,TC*U0/L 
3      UOO  =  UNN 

DEOSOO  =  DEOSNN 
2    CONTINUE 
1  CONTINUE 
C 
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2001  FORMAT  (215) 
END 


SUBROUTINE  DSTABL(DEOS,WR,WI,OMEGA) 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
DIMENSION  WR(4),WI(4) 
DEOS=-1.0D+20 
DO  1  1=  1,4 

IF  (WR(I).LT.DEOS)  GO  TO  1 

DEOS  =  WR(I) 

IJ  =  I 
1  CONTINUE 
OMEGA  =  WI(IJ) 
OMEGA  =  DABS(OMEGA) 
RETURN 
END 

SUBROUTINE  TRAP(N,A,B,OUT) 

NUMERICAL  INTEGRATION  ROUTINE  USING  THE  TRAPEZOIDAL  RULE 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
DIMENSION  A(1),B(1) 

Nl  =N-1 
OUT  =  0.0 
DO  1  I  =  1,N1 

OUT1  =  0.5*(A(I)+A(I+1))*(B(I+1)-B(I)) 

OUT  =  OUT+OUTl 
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1    CONTINUE 
RETURN 
END 

C     FUNCTIONS  USED  IN  NEWTON  RAPHSON  ROUTINE 
FUNCTION  FUNC(THETAO) 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
COMMON  P1,P2,P3,P4 

FUNC=  P1*DTAN(THETA0)  +  P2*DCOS(THETA0)  +  P3*DSIN(THETA0) 
&       +  P4*(DTAN(THETA0))**2 
RETURN 
END 

FUNCTION  DFUNC(THETAO) 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
COMMON  PI, P2,P3,P4 

DFUNC=P1*(1.DO/DCOS(THETAO))**2  -  P2*DSIN(THETA0)  + 
&    P3*DCOS(THETA0)  +  P4*2.D0*DTAN(THETA0)*(l.D0/DCOS(THETA0))**2.D0 
RETURN 
END 
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